set.seed(123)
d <- sample(100, 10)
d [1] 31 79 51 14 67 42 50 43 97 25
Omar E. Barrantes Sotela
14 de septiembre de 2022
La autocorrelación espacial es un concepto importante en las estadísticas espaciales. Su fundamento teórico es lo que permite la interpolación espacial. En muchas ocasiones su calculo y propiedades son malentendidas.
La autocorrelación (ya sea espacial o no) es una medida de similitud (correlación) entre observaciones cercanas. Para comprender la autocorrelación espacial, primero es útil considerar la autocorrelación temporal.
Si mide algo sobre el mismo objeto a lo largo del tiempo, por ejemplo, el peso de una persona, es probable que dos observaciones que estén cercanas entre sí en el tiempo también sean similares en la medición. Digamos que en un par de años su peso pasó de \(50\) a \(80\) kg. Es poco probable que fuera \(60\) kg un día, \(50\) kg el siguiente y \(80\) al día siguiente. Por el contrario, probablemente subió gradualmente, con la reducción gradual ocasional, o incluso revertir en la dirección. Lo mismo puede ser cierto con su cuenta bancaria, pero eso también puede tener una tendencia mensual marcada. Para medir el grado de asociación a lo largo del tiempo, podemos calcular la correlación de cada observación con la siguiente observación.
En este caso d es un vector de observaciones diarias.
La autocorrelación calculada anteriormente es muy pequeña. Aunque esta es una muestra aleatoria, (casi) nunca obtiene un valor de cero. Calculamos la autocorrelación de “un retraso”, es decir, comparamos cada valor con su vecino inmediato y no con otros valores cercanos.
Después de ordenar los números en la variable d, la autocorrelación se vuelve muy fuerte (como era de esperar).
La función acf muestra la autocorrelación calculada de una manera ligeramente diferente para varios retrasos (es 1 para cada punto en sí mismo, muy alta cuando se compara con el vecino más cercano, y se reduce gradualmente).
El concepto de autocorrelación espacial es una extensión de la autocorrelación temporal. Aunque es un poco más complicado. El tiempo es unidimensional, y solo va en una dirección, siempre adelante. Los objetos espaciales tienen (al menos) dos dimensiones y formas complejas, y puede no ser obvio cómo determinar qué está “cerca”.
Las medidas de autocorrelación espacial describen el grado en que las observaciones (valores) en ubicaciones espaciales (ya sean puntos, áreas o celdas de trama) son similares entre sí. Entonces necesitamos dos cosas: observaciones y ubicaciones.
La autocorrelación espacial en una variable puede ser exógena (es causada por otra variable espacialmente autocorrelacionada, por ejemplo, lluvia) o endógena (es causada por el proceso en juego, por ejemplo, la diseminación de una enfermedad).
Una estadística comúnmente utilizada que describe la autocorrelación espacial es la I de Moran, y discutiremos esto aquí en detalle. Otros índices incluyen C de Geary y, para datos binarios, el índice de recuento de uniones. El semi-variograma también expresa la cantidad de autocorrelación espacial en un conjunto de datos (el cual se estudiará en el tema de interpolación) (Barrantes-Sotela & Solano-Mayorga, 2020).
library(rgdal)
library(raster)
#Carga un shapefile en R
shape <- "./shp/San_Pedro_CU.shp"
p <- readOGR(dsn = shape,layer = 'San_Pedro_CU')OGR data source with driver: ESRI Shapefile
Source: "G:\UNA\Clases\ECG-UNA\MC\gei411-mc\demo\demo03\shp\San_Pedro_CU.shp", layer: "San_Pedro_CU"
with 14 features
It has 28 fields
Integer64 fields read as strings: ugm tot_viv tot_hog hog_nbi tot_pp pp_noCCSS
[1] "+proj=tmerc +lat_0=0 +lon_0=-84 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs"
Sí se supone el interés en la autocorrelación espacial en la variable “Necesidades Básicas Insatisfechas”. Si hubiera autocorrelación espacial, las Unidades Mínimas Geocensalescon valor similar estarían agrupadas espacialmente.
Se construye un diagrama de los polígonos. Se utiliza la función coordinates para obtener los centroides de los polígonos para colocar las etiquetas.
Ahora necesitamos determinar qué polígonos están “cerca” y cómo cuantificar eso. Aquí usaremos la adyacencia como criterio. Para encontrar polígonos adyacentes, podemos usar el paquete spdep.
[1] "nb"
Neighbour list object:
Number of regions: 14
Number of nonzero links: 64
Percentage nonzero weights: 32.65306
Average number of links: 4.571429
Link number distribution:
3 4 5 6 8
4 4 2 3 1
4 least connected regions:
O M A B with 3 links
1 most connected region:
K with 8 links
summary(w) genera un resumen de los vecinos. El número promedio de vecinos (polígonos adyacentes) es 4.57, 4 con 3 enlaces, 4 con 4 enlaces, 2 con 5 enlaces, 3 con 6 enlaces y 1 con 8 (¿Cúal es ese? Respuesta: k.).
Para mayores detalles se puede observar la estructura (str) de w.
List of 14
$ : int [1:6] 2 3 7 8 9 14
$ : int [1:4] 1 7 8 14
$ : int [1:4] 1 8 9 14
$ : int [1:3] 5 7 8
$ : int [1:5] 4 6 7 8 9
$ : int [1:3] 5 8 9
$ : int [1:5] 1 2 4 5 8
$ : int [1:8] 1 2 3 4 5 6 7 9
$ : int [1:6] 1 3 5 6 8 14
$ : int [1:4] 11 12 13 14
$ : int [1:4] 10 12 13 14
$ : int [1:3] 10 11 13
$ : int [1:3] 10 11 12
$ : int [1:6] 1 2 3 9 10 11
- attr(*, "class")= chr "nb"
- attr(*, "region.id")= chr [1:14] "H" "I" "G" "O" ...
- attr(*, "call")= language poly2nb(pl = p, row.names = p$cod)
- attr(*, "type")= chr "queen"
- attr(*, "sym")= logi TRUE
integer(0)
Podemos transformar w en una matriz de ponderaciones espaciales. Una matriz de ponderaciones espaciales refleja la intensidad de la relación geográfica entre las observaciones.
| H | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 |
| I | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 |
| G | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 |
| O | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| N | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 |
| M | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 |
| L | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| K | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| J | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 |
| D | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
| C | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 |
| A | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 0 |
| B | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 |
| F | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 |
El cálculo de la I’ Moran corresponde a la siguiente formula
\[ I=\frac{n} {\sum_{i=1}^{n} {(y_i - \bar{y})^2}} \frac{\sum_{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}})} {\sum_{i}\sum_{j}w_{ij}} \]
Otra forma de verlo sería
\[ I={\frac {n}{w}}{\frac {\sum _{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}})}{\sum _{i}(y_{i}-{\bar {y}})^{2}}} \]
I’ Moran es una versión expandida de la fórmula para calcular el coeficiente de correlación. Lo principal que se agregó es la matriz de ponderaciones espaciales.
Primero se determina el número de observaciones
Se obtiene los valores \(y\) y \(\bar{y}\)
Ahora se obtiene \((y_i - \bar{y})(y_j - \bar{y})\), para todos los pares de datos. Se puede realizar con dos métodos.
Se procede a obtener una matriz de los pares multiplicados
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 24.566184 -24.799844 30.018256 -39.322180 -62.914780 72.247027
[2,] -24.799844 25.035727 -30.303773 39.696191 63.513191 -72.934202
[3,] 30.018256 -30.303773 36.680327 -48.049109 -76.877709 88.281098
[4,] -39.322180 39.696191 -48.049109 62.941556 100.705356 -115.643137
[5,] -62.914780 63.513191 -76.877709 100.705356 161.126756 -185.026937
[6,] 72.247027 -72.934202 88.281098 -115.643137 -185.026937 212.472270
[7,] -43.089066 43.498906 -52.651994 68.971070 110.352470 -126.721223
[8,] -62.914780 63.513191 -76.877709 100.705356 161.126756 -185.026937
[9,] -62.914780 63.513191 -76.877709 100.705356 161.126756 -185.026937
[10,] 102.282984 -103.255844 124.983056 -163.720580 -261.949980 300.805427
[11,] 43.301484 -43.713344 52.911556 -69.311080 -110.896480 127.345927
[12,] 24.566184 -24.799844 30.018256 -39.322180 -62.914780 72.247027
[13,] 1.023148 -1.032880 1.250220 -1.637716 -2.620316 3.008991
[14,] -2.049837 2.069334 -2.504766 3.281098 5.249698 -6.028394
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] -43.089066 -62.914780 -62.914780 102.282984 43.301484 24.566184
[2,] 43.498906 63.513191 63.513191 -103.255844 -43.713344 -24.799844
[3,] -52.651994 -76.877709 -76.877709 124.983056 52.911556 30.018256
[4,] 68.971070 100.705356 100.705356 -163.720580 -69.311080 -39.322180
[5,] 110.352470 161.126756 161.126756 -261.949980 -110.896480 -62.914780
[6,] -126.721223 -185.026937 -185.026937 300.805427 127.345927 72.247027
[7,] 75.578184 110.352470 110.352470 -179.404266 -75.950766 -43.089066
[8,] 110.352470 161.126756 161.126756 -261.949980 -110.896480 -62.914780
[9,] 110.352470 161.126756 161.126756 -261.949980 -110.896480 -62.914780
[10,] -179.404266 -261.949980 -261.949980 425.862184 180.288684 102.282984
[11,] -75.950766 -110.896480 -110.896480 180.288684 76.325184 43.301484
[12,] -43.089066 -62.914780 -62.914780 102.282984 43.301484 24.566184
[13,] -1.794602 -2.620316 -2.620316 4.259948 1.803448 1.023148
[14,] 3.595413 5.249698 5.249698 -8.534637 -3.613137 -2.049837
[,13] [,14]
[1,] 1.02314847 -2.04983724
[2,] -1.03288010 2.06933418
[3,] 1.25021990 -2.50476582
[4,] -1.63771582 3.28109847
[5,] -2.62031582 5.24969847
[6,] 3.00899133 -6.02839439
[7,] -1.79460153 3.59541276
[8,] -2.62031582 5.24969847
[9,] -2.62031582 5.24969847
[10,] 4.25994847 -8.53463724
[11,] 1.80344847 -3.61313724
[12,] 1.02314847 -2.04983724
[13,] 0.04261276 -0.08537296
[14,] -0.08537296 0.17104133
Ahora es posible multiplicar esta matriz con los pesos para establecer en cero el valor para los pares que no son adyacentes.
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
H 0.000000 -24.799844 30.018256 0.00000 0.0000 0.0000 -43.08907
I -24.799844 0.000000 0.000000 0.00000 0.0000 0.0000 43.49891
G 30.018256 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
O 0.000000 0.000000 0.000000 0.00000 100.7054 0.0000 68.97107
N 0.000000 0.000000 0.000000 100.70536 0.0000 -185.0269 110.35247
M 0.000000 0.000000 0.000000 0.00000 -185.0269 0.0000 0.00000
L -43.089066 43.498906 0.000000 68.97107 110.3525 0.0000 0.00000
K -62.914780 63.513191 -76.877709 100.70536 161.1268 -185.0269 110.35247
J -62.914780 0.000000 -76.877709 0.00000 161.1268 -185.0269 0.00000
D 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
C 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
A 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
B 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
F -2.049837 2.069334 -2.504766 0.00000 0.0000 0.0000 0.00000
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
H -62.91478 -62.914780 0.000000 0.000000 0.000000 0.000000 -2.049837
I 63.51319 0.000000 0.000000 0.000000 0.000000 0.000000 2.069334
G -76.87771 -76.877709 0.000000 0.000000 0.000000 0.000000 -2.504766
O 100.70536 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
N 161.12676 161.126756 0.000000 0.000000 0.000000 0.000000 0.000000
M -185.02694 -185.026937 0.000000 0.000000 0.000000 0.000000 0.000000
L 110.35247 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
K 0.00000 161.126756 0.000000 0.000000 0.000000 0.000000 0.000000
J 161.12676 0.000000 0.000000 0.000000 0.000000 0.000000 5.249698
D 0.00000 0.000000 0.000000 180.288684 102.282984 4.259948 -8.534637
C 0.00000 0.000000 180.288684 0.000000 43.301484 1.803448 -3.613137
A 0.00000 0.000000 102.282984 43.301484 0.000000 1.023148 0.000000
B 0.00000 0.000000 4.259948 1.803448 1.023148 0.000000 0.000000
F 0.00000 5.249698 -8.534637 -3.613137 0.000000 0.000000 0.000000
attr(,"call")
nb2mat(neighbours = w, style = "B")
Ahora se suman los valores, para obtener un aproximado del I de Moran:
\[ \sum _{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}}) \]
El siguiente paso es dividir este valor por la suma de los pesos.
Se calcula la varianza inversa de y
El último paso para calcular el I de Moran.
Esta es una forma simple (pero cruda) de estimar el valor esperado de I. de Moran, es decir, el valor que obtendría en ausencia de autocorrelación espacial (si los datos fueran espacialmente aleatorios). Por supuesto, nunca se espera eso, pero así es como se hace en las estadísticas. Se debe tener en cuenta que el valor esperado se acerca a cero si n se vuelve grande, pero que no lo suficiente para valores pequeños de n.
Después de hacer esto a mano, ahora se usará el paquete spdep para calcular el I de Moran y realizar una prueba de significancia. Para hacer esto, se necesita crear un objeto de pesos espaciales tipo listw (en lugar de la matriz que utilizamos anteriormente). Para obtener el mismo valor que el anterior, usamos style = 'B' para usar pesos de distancia binarios (TRUE / FALSE).
Characteristics of weights list object:
Neighbour list object:
Number of regions: 14
Number of nonzero links: 64
Percentage nonzero weights: 32.65306
Average number of links: 4.571429
Weights style: B
Weights constants summary:
n nn S0 S1 S2
B 14 196 64 128 1288
Ahora podemos usar la función moran. Echa un vistazo con ?moran. La función se define como moran (y, ww, n, Szero (ww)). Tenga en cuenta los argumentos impares n y S0. Creo que son raros, porque “ww” tiene esa información. De todos modos, los suministramos y funciona. Probablemente haya casos en los que tenga sentido usar otros valores.
$I
[1] 0.1609378
$K
[1] 2.158294
[1] 64
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
H 0.000000 -24.799844 30.018256 0.00000 0.0000 0.0000 -43.08907
I -24.799844 0.000000 0.000000 0.00000 0.0000 0.0000 43.49891
G 30.018256 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
O 0.000000 0.000000 0.000000 0.00000 100.7054 0.0000 68.97107
N 0.000000 0.000000 0.000000 100.70536 0.0000 -185.0269 110.35247
M 0.000000 0.000000 0.000000 0.00000 -185.0269 0.0000 0.00000
L -43.089066 43.498906 0.000000 68.97107 110.3525 0.0000 0.00000
K -62.914780 63.513191 -76.877709 100.70536 161.1268 -185.0269 110.35247
J -62.914780 0.000000 -76.877709 0.00000 161.1268 -185.0269 0.00000
D 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
C 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
A 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
B 0.000000 0.000000 0.000000 0.00000 0.0000 0.0000 0.00000
F -2.049837 2.069334 -2.504766 0.00000 0.0000 0.0000 0.00000
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
H -62.91478 -62.914780 0.000000 0.000000 0.000000 0.000000 -2.049837
I 63.51319 0.000000 0.000000 0.000000 0.000000 0.000000 2.069334
G -76.87771 -76.877709 0.000000 0.000000 0.000000 0.000000 -2.504766
O 100.70536 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
N 161.12676 161.126756 0.000000 0.000000 0.000000 0.000000 0.000000
M -185.02694 -185.026937 0.000000 0.000000 0.000000 0.000000 0.000000
L 110.35247 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
K 0.00000 161.126756 0.000000 0.000000 0.000000 0.000000 0.000000
J 161.12676 0.000000 0.000000 0.000000 0.000000 0.000000 5.249698
D 0.00000 0.000000 0.000000 180.288684 102.282984 4.259948 -8.534637
C 0.00000 0.000000 180.288684 0.000000 43.301484 1.803448 -3.613137
A 0.00000 0.000000 102.282984 43.301484 0.000000 1.023148 0.000000
B 0.00000 0.000000 4.259948 1.803448 1.023148 0.000000 0.000000
F 0.00000 5.249698 -8.534637 -3.613137 0.000000 0.000000 0.000000
attr(,"call")
nb2mat(neighbours = w, style = "B")
[1] 132
Ahora podemos probar la significancia. Primero analíticamente, usando lógica basada en regresión lineal y suposiciones.
Moran I test under normality
data: p$hog_NBI_p
weights: ww
Moran I statistic standard deviate = 1.7582, p-value = 0.03935
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.16093780 -0.07692308 0.01830159
En lugar del enfoque anterior, debe usar la simulación de Monte Carlo. Ese es el método preferido (de hecho, el único método). Funciona de manera que los valores se asignan aleatoriamente a los polígonos, y se calcula el I de Moran. Esto se repite varias veces para establecer una distribución de valores esperados. El valor observado de I de Moran se compara con la distribución simulada para ver qué tan probable es que los valores observados puedan considerarse aleatorios.
Se puede hacer un “diagrama de dispersión de Moran” para visualizar la autocorrelación espacial. Primero obtenemos los valores vecinos para cada valor.
Se remueven los valor iguales a 0.
Se calculan los valores de los promedios de los vecinos.
ams <- aggregate(ms[,2:3], list(ms[,1]), FUN=mean)
ams <- ams[,-1]
colnames(ams) <- c('y', 'y espacialmente retrasado')
head(ams) y y espacialmente retrasado
1 17.65 10.68000
2 7.69 11.31000
3 18.75 14.96500
4 4.76 4.00000
5 0.00 12.01000
6 4.00 10.03333
Para finalmente graficar
plot(ams)
reg <- lm(ams[,2] ~ ams[,1])
abline(reg, lwd=2)
abline(h=mean(ams[,2]), lt=2)
abline(v=ybar, lt=2)Tenga en cuenta que la pendiente de la línea de regresión:
es casi el mismo valor que el I. de Moran
Aquí hay un enfoque más directo para lograr lo mismo (pero es de esperar que lo anterior aclare cómo se calcula esto en realidad). Tenga en cuenta la estandarización de filas de la matriz de ponderaciones:
rwm <- mat2listw(wm, style='W')
# Revisa si las filas suman hasta 1
mat <- listw2mat(rwm)
apply(mat, 1, sum)[1:15] H I G O N M L K J D C A B F <NA>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 NA
Se vuelve a graficar
Analice el gráfico correspondiente
¿Cuál es su opinión con respecto a este método?
¿Tiene importancia el análisis espacial de fenómenos mediante el I’ Moran?